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In this paper 1 describe a simple numerical procedure to compute synthetic horizon altitude profiles for any given site. 
The method makes use of a simplified model of local Earth's curvature, and it is based on the availability of digital 
elevation models describing the topography of the area surrounding the site under study. Examples constructed using the 
Shuttle Radar Topographic Mission (SRTM) data (with 90m horizontal resolution) are illustrated, and compared to direct 
theodolite measurements. The proposed method appears to be reliable and applicable in all cases when the distance to the 
local horizon is larger than ~10 km, yielding a rms accuracy of ~0. 1 degrees (both in azimuth and elevation). Higher 
accuracies can be achieved with higher resolution digital elevation models, like those produced by many modern national 
geodetic surveys. 
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1 Introduction 

When studying the orientation of buildings, tombs, or any 
other man made structure, having an handy description of 
the natural horizon is a fundamental step. In most of the 
cases this is done directly measuring the horizon altitude at 
the relevant azimuths using, for instance, a theodolite and 
sun fixes. Depending on the angular sampling one wants to 
achieve, this might turn into a rather long and boring pro- 
cedure. An alternative possibility is the calibration of a few 
points using direct theodolite readings coupled to digital, 
rectified photography. However, this might turn to be diffi- 
cult, or even impossible, in the case the natural horizon is 
not visible because of modern constructions or vegetation. 

In this paper I propose an alternative solution, which is 
based on a simple geometrical model for Earth's local shape 
and the availability of a digital elevation model (DEM) for 
the area under examination. The idea is rather simple. For a 
given observer's site, the line of site (hereafter LOS) eleva- 
tion profile (LOSEP) along the input azimuth is extracted. 
Then, for each point along the LOS, the apparent altitude 
above the local horizontal is estimated and the maximum 
is found. By definition, this is the natural horizon altitude 
as seen from the observing site. Repeating the same proce- 
dure for all azimuths (with a given angular step) will finally 
allow one to retrieve the horizon altitude profile. The angu- 
lar resolution clearly depends on the horizontal sampling in 
the DEM, while the accuracy is related to the DEM vertical 
accuracy. 

Since the available DEMs provide the elevation above 
a given ellipsoid (typically the WGS84), the first thing one 
needs to take into account is Earth's curvature, which makes 



a far mountain appear lower than it would if the Earth's sur- 
face were flat. Then, although this is a second order effect, 
one has to coiTect for the teiTestiial refraction (which has 
the opposite effect). These two coiTections are discussed in 
the next two sections. 

2 Correcting for Earth's curvature 

Let us consider an observer placed in A at an elevation Ka 
(above the sea level), and a point B (elevation hs), located 
at a distance d from A. If Earth were flat, the altitude a of 
B seen from A would simply be: 

(hs-hA^ 

a — arctan ; — 

V d 

However, because of Earth's curvature, far objects sub- 
tend altitudes which are smaller than those given by the pre- 
vious expression. Using a simplified model for Earth's local 
curvature, one can derive the following approximate expres- 
sion (see Appendix [a): 

hs — hA Id 
d 2~R 



a ~ arctan 



(1) 
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where R is the local radius of curvature (see Equation Al i. 
For example, suppose that a /ib=3.0 km mountain peak is 
observed from a site placed d=\QQ km away, at an elevation 
/iA=0.1 km. Using i?=6370 km one obtains q;=1°.21. Ne- 
glecting Earth's curvature one would get a=l°.66, almost 
half a degree off. It can be easily shown that for ellipsoidal 
distances d <300 km the error introduced by the approxi- 
mated expression (fTl) is less than 0.01 degrees, which is well 
below the typical archaeoastronomical needs. 

In the cases where \hB — hA\/d <^ 1, using Equationfl] 
one can show that the apparent altitude decreases at a rate 
of 1/2R, which is about 16".2 km^i (or 0°.0045 km^^). 
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Fig. 1 Example LOSEP extraction for the Roman town 

Aquileia (A=13°.367 E, (/)= 45°.767 N), with ^=339°. The 
underlying DEM is from the SRTM90 data set, and is 
displayed in the Universal Transverse Mercator projection 
(UTM, 33N). The white area in the lower part of the map is 
the Adriatic Sea. 

3 Correcting for terrestrial refraction 

As the light travels across Earth's atmosphere it bends, in 
such a way that a far object appears to be higher than actu- 
ally is. This phenomenon is known as terrestrial refractiorT] 
If a is the unrefracted horizon altitude (corrected for 
Earth's curvature), then the apparent horizon altitude a' is 
simply given by a' = a + Rt- Several analytical descrip- 
tions of Rt have been proposed, but the one developed by 



Bomsford ( 19801 appears to be the most accurate (see for 
instance Sampson et al. |2003[ ). In this formulation the ter- 
restrial refraction is given by 

Rt = Kd , 

where d is the distance between the observer an the natural 
horizon (in km), and k (degrees km~^) is defined as: 

1 0.252P / ^ „ dT\ 180 



(2) 



where R is Earth's radius (in km), P is the atmospheric 
pressure (in millibars), T is the ground air temperature (in 
K), and dT/dz is the ground atmospheric vertical temper- 
ature gradient (in K km^^). For typical atmospheric condi- 
tions (P=1000 mb, T=293 K, dT/dz=-\0 K km'^) one 
finds that k ^2.3 arcsec km~^ (or 0.00064 degrees km^^). 
This is about seven times smaller than the effect produced 



' As opposed to the so called astronomical refraction. The phenomenon 
is very similar, but in that case the light rays have to cross the whole atmo- 
sphere, which makes its description much more complex, especially when 
one is to consider objects very close to the horizon. 



by Earth's curvature. In the case of the example discussed 
in the previous section {d=\QQ km) this would turn into a 
correction of ^0.06 degrees. 

One can simplify the previous formula introducing the 
refraction constant K: 



K = 0.273 34.2 



dT 

dz 

With this setting, k (in degrees km^^) can be written as: 



K = 0.00829 K 



P 



(3) 



Typical values of K range from ^5 (during the day) 
to ^10 (at sunset/sunrise or night). The fluctuation of at- 
mospheric conditions introduce a variation in k, which can 
range from ^1.5 to ^4 arcsec km^^. For instance, on a 
distance of 100 km this translates into a peak-to-peak ho- 
rizon altitude change of ^0.07 degrees. Therefore, even in 
the hypothesis all other quantities are known with negligi- 
ble errors, this sets the natural maximum accuracy one can 
achieve on the horizon apparent altitude. 

4 Digital Elevation Model 

The proposed method is based on the availability of a DEM, 
i.e. an array of values giving the elevation of a given point 
above the underlying ellipsoid. If A and are the longitude 
and latitude of a point P on the ellipsoid, I will express its 
elevation z as z = ^(A, 0) or, alternatively, as z = g{P)- 

Because of its nature, a DEM is a discrete collection 
of data points, obtained with some spatial resolution. Typi- 
cally, the data are distributed on a regular grid, with a hori- 
zontal sampling that 1 will indicate as AL This corresponds 
to the minimum scale one can resolve in the DEM. Obvi- 
ously, if one is to get the elevation of a point which does 
not coincide with one of the grid nodes, then one will have 
to use some interpolation method (nearest neighbor, linear 
interpolation, bi-cubic spline interpolations, etc). No matter 
what method is used, though, the resolution is dictated by 
the sampling. In the following I will assume that the DEM 
is given with a regular sampling AZ, which is to say that 
the data points have been obtained/re-gridded at a constant 
step. Also, I will indicate with a^ the rms uncertainty on 
each DEM data point. 

Each country has its own geographical survey program- 
me, aiming at mapping its territory with a certain resolu- 
tion. These data may or may not be available to the reader. 
Therefore, in this article 1 will consider the Shuttle Radar 
Topographic Mission, which has the advantage of having a 
relatively good horizontal resolution, a good vertical accu- 
racy, an almost world wide coverage, and, most importantly, 
is freely availablq^ 



The data can be freely downloaded at the following URL: 



|http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp| 
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4.1 The Shuttle Radar Topographic Mission 

The Shuttle Radar Topographic Mission (SRTM) provides 
an homogeneous coverage of Earth's elevation (Farr et al. 
2007[ ). For the US territory it is distributed with a hori- 
zontal resolution of 30 m (SRTM30), while for the rest of 
the planet data were averaged within bins of 90x90 m^ 
(SRTM90), hence including 3x3 original data points each. 
The 90%-level absolute vertical accuracy of the SRTM orig- 
inal data is better than 9 m (Farr et al. |2007| l, which corre- 
sponds to an rms accuracy of 5.5 m. Therefore, the formal 
rms error on the 90 m data is expected to be a^ ^1.8 m. I 
note, however, that this is strictly true only if the terrain is 
smooth on scales smaller than 30 m, otherwise the averag- 
ing process within 90 m bins produces an artificial smooth- 
ing, turning into much larger deviations from the real values. 
Clearly this is going to affect very steep mountain regions, 
where the elevation can change significantly on scales of 
tens of meters. As a consequence, sharp mountain peaks ap- 
pear in the SRTM90 data with lower than real elevations. 

The altitude uncertainty Ua (in degrees) implied by an 
elevation uncertainty g^ can be estimated through the fol- 
lowing approximate expression: 

180%/2 a^ 

o-Q ~ -r , 

TT a 

where d is the distance between the observer and the point 
being surveyed, expressed in the same units as a 2- For in- 
stance, the formal elevation error of the SRTM90 data (1.8 
m) translates into an altitude error of about 0. 1 degrees at 
d ~1.5 km. 

As for the SRTM horizontal accuracy, this is better than 
20 m (at the 90%-level; Farr et al. 2007| l. This corresponds 
to an rms deviation (T(=14 m). Therefore, in the SRTM90 
case, this is much better than the horizontal resolution Al. 
As a consequence, the rms error on the azimuth 1^9 of a given 
data point is approximately 



T, 



180 C7i 



TT 



kn 



At (i=10 km this corresponds to an uncertainty of ~0.1 
degrees, which is a fifth of the apparent diameter of the sun 
(and the moon). Obviously, if the horizon is farther away, 
the azimuth uncertainty becomes proportionally smaller If 
one is after higher azimuth accuracies, or the horizon in the 
direction of interest is closer than 10 km, then a DEM with 
better sampling is requireq^ 

In this respect, I note that a prehminary version of the 
data produced by the Advanced Spaceborne Thermal Emis- 
sion and Reflection Radiometer (ASTERrll was recently re- 
leased. Nominally they have an horizontal resolution of 30 
m, which is a factor 3 better than the SRTM90. However, a 



' In many countries surveys with horizontal resolutions of 15 m or better 
are either available or under construction. The reader should get in touch 
with the local authoiities to verify the public availability of higher resolu- 
tion data. 

* See: http://asterweb.jpl.nasa.gov/index.asp 
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Fig. 2 Upper panel: LOSEP along ip=339° for the Roman 
town of Aquileia (see Figure [T]). Lower panel: the corre- 
sponding altitude profile. The blue arrow marks the max- 
imum altitude, corresponding to the natural horizon along 
the chosen direction. 

close look to the data shows that the effective horizontal res- 
olution is actually close to 90 m and the vertical accuracy is 
not as good as the one of SRTM, at least not in the released 
version. 

5 LOSEP extraction and horizon 
determination 

Having a DEM at hand, one can extract the LOSEP start- 
ing from the observer's position Pq along a given azimuth 
V? as a function of distance d from the observer. This corre- 
sponds to solving the forward geodesic problem, i.e. com- 
puting the running end point P of a geodesic path on the 
ellipsoid, given the start point Pq, a path length d, and a 
starting azimuth (p. As is standard in geodesy, this is done 
numerically, using the iterative algorithm devised by Vin- 



centy ( 1975 1. Although the solution is non-analytical, for 
the sake of clarity I will indicate it as 

P^f{Po,d,ip). 

Given the discrete nature of the DEM, the LOSEP is ex- 
tracted at a number of points, which are separated by some 
constant length. Obviously, the maximum resolution is at- 
tained when this length is equal to the DEM resolution AL 

With these definitions, the LOSEP extraction is done 
through the following steps: 

1. elevation zq = g{Po) is extracted from the DEM at 

(Ao,(/'o); 

2. the path length for the i-th step along the LOS is updated 

to di = Al i; 
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Fig. 3 Comparison between the theodolite readings (empty circles) and the synthetic profiles (W to the left and E to the 
right) computed using SRTM (red) and ASTER (blue) DEMs for the site of St. Martin (Artegna - Italy). Top panel: horizon 
distance. Middle panel: horizon profiles. Bottom panel: Aa residuals (computed— measured). The errorbars indicate the 
5-cr error of the synthetic profile. 



3. the coordinates (Ai, 4>i) of the running point P are com- 
puted as /(Pq, d^, (^) 

4. the elevation of the running point is computed as zi = 

5. the apparent altitude ai of P as seen from Pq is cal- 
culated with Equation |A4| and corrected for terrestrial 
refraction; 

6. the cycle is repeated from stepl2j until di reaches a max- 
imum value dmax- 

Once this is done, the horizon elevation along azimuth ip 
is determined as a{if) — maxjofi}. I note that in most of the 
cases it is sufficient to consider distances within dmax —200 
km. The exact value depends on the maximum elevation in 
the area surrounding the site under study. From Equation [T] 
one has that the distance at which a mountain peak appears 
at an altitude a=0 is given by d w ^2RAz, where Az is the 
difference in elevation. For Az <4.0 km (which is true for 
most of the sites on the planet), it is d <225 km. Inciden- 
tally, here emerges one of the advantages of horizon syn- 
thesis over direct measurements. If the horizon is far away 
(d >50-100 km), it might be very difficult, if not impossi- 
ble, to have a sufficient atmospheric transparency to be able 
to actually see it. 

An example LOSEP extraction is presented in Figure [T] 
where I have chosen the Roman town Aquileia (Italy; fo- 
unded in 181 BCE) as the observing site. For the sake of the 
example I have traced the LOS along one of the two cardi- 
nal directions of the town [ip ~339°), in the NW direction. 
The site is marked with a circle, while the selected line of 
sight is traced by a solid line. The underlying DEM (from 
the SRTM90 data set) has been projected to the Universal 



Transverse Mercator system (see Hager, Behensky & Drew 
[1989). 

The LOSEP extracted along this direction is shown in 
Figure |2] (upper panel), and it reaches about 3000 m in the 
Austrian Alps, at a distance of 140 km from Aquileia. How- 
ever, the horizon is located at about 100 km, at an elevation 
of ^2600 m, and it subtends an apparent altitude of ^l°.l 
(Figurel2] lower panel). 

Clearly the whole process can be repeated for (p ranging 
from 0° to 360° with a given step Aip. This will finally give 
the full horizon profile a{(p). 

6 Validation 

The simplest way of validating the procedure outlined in the 
previous section is a comparison between the synthetic pro- 
file and direct horizon measurements. The natural horizon 
altitude can be determined using a theodolite and sun fixes, 
taking measurements with some azimuth step. The typical 
accuracy one can achieve with this method is of the order of 
a tenth of a degree or better, depending on the quality of the 
instrument useq3 

In the course of studying the orientation of the ancient 
church of St. Martin (Artegna - Italy; 13°. 1528 E, 46°.2415 
N, 267 m a.s.l.), I had measured the natural horizon alti- 
tude in a range around E and W directions. The building 
is located on top of St. Martin Hill, about 50 m above the 



^ It must be noticed that when the horizon is close, trees can substan- 
tially affect direct measurements, producing systematically higher theodo- 
lite readings. 
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surrounding plain. The E horizon is dominated by the pres- 
ence of Mt. Faet (750 m a.s.L), whose top is located at less 
than 2 km from St. Martin. Seen from the site, the top sub- 
tends an angle of about 14 degrees. On the contrary, the W 
horizon is located at more than 10 km, and has an altitude 
around 3 degrees. The different distances to the E and W 
horizon locations from the observing site make this an ideal 
test case. The comparison between the theodolite readings 
and the synthetic profile is presented in Fig. [3] 

In the case of the W horizon, the deviations are always 
smaller than 0.2 degrees, and in most of the cases they are 
less than 0.1 degrees, i.e. well within the expected rms er- 
rors. In general, the SRTM DEM gives a better result than 
the ASTER one (Fig. [3] left panels). Things are different 
for the E horizon, for which the match is good only along 
the smooth declining ridge of Mt. Faet in the azimuth range 
68°-78°, where the horizon is at a distance of about 1.9 km 
(Fig. [31 right panels). Between 64° and 68°, where the dis- 
tance drops to about 1.2 km, the deviations exceed 0°.5. But 
the largest discrepancies are seen for azimuths larger than 
84°, where they are larger than 1°.5. This very clearly illus- 
trates the effect of a close horizon coupled to a steep slope, 
where the smoothing intrinsic to the DEM data produces 
systematically lower elevations and, in turn, altitudes. 

These two examples demonstrate that the procedure de- 
scribed here, when used in conjunction to SRTM90 or AS- 
TER DEMs, gives results accurate to ^0.1 degrees only 
for horizon distances larger than ~10 km. This does not 
exclude that acceptable results may be achieved also for 
shorter distances, but this will depend on the horizon mor- 
phology, and cannot be assumed a priori. 



7 Example application: the case of Maes 
Howe 

The method described in this paper can be applied any time 
there is the need of having a handy description of the nat- 
ural horizon in an archaeo-astronomical analysis. To show 
the potential of the technique, I close this article with an 
application to a famous and well studied case, where the 
horizon shape is playing a fundamental role: Maes Howe. 

Maes Howe is a neolithic site located on Mainland, Ork- 
ney, Scotland (3°. 1879 W, 58°.9981 N, 20 m a.s.L). Besides 
being the site of a chambered cairn and passage grave, Maes 
Howe offers a spectacular event. Around twenty days before 
(and after) the winter solstice, the sun sets below the horizon 
defined by Ward Hill, placed at about 14 km SW of Maes 
Howe. However, after 7-8 minutes, it reappears for a couple 



of minutes, just before definitely setting (Reijs 1998 2001 
See also Magli ,2009, his Figs. 3 and 4). This is due to a 
combination of the sun apparent path and the morphology 
of the horizon. 

This phenomenon is very well reproduced using the syn- 
thetic horizon computed for Maes Howe, once one applies 
a standard astronomical refraction correction to the posi- 
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Fig. 4 Lower panel: horizon location for Maes Howe. The 
map is a contour plot of the SRTM DEM. The position of 
Ward Hill is marked by a blue circle. Upper panel: appar- 
ent path of the sun on December 1st (Sq=—21° .9; solid 
red line). The dashed line traces the real path of the sun 
(see text). The red circles mark the position of the sun (to 
scale) with a time interval of 2 minutes. The blue curve is 
the SRTM synthetic horizon as seen from Maes Howe. 



tion of the surj^ This is illustrated in Fig. |4] where I have 
traced the apparent trajectory of the sun over the horizon. 
The sun reappears for a couple of minutes at an azimuth 
around 223°, in full agreement with direct measurements 
(Reijs [T998) , when the sun declination is 5q=— 21° .9 (cor- 
responding to December 1st of today's calendar). An ac- 
curate reproduction is obtained also for the sun reappea- 
rance in the direction of Kame of Hoy, which takes place for 
Sq — — 17°. 3 (February 1st), at an azimuth of about235°. 5. 
Since in 2700 BCE (the supposed date for the construction 
of Maes Howe) the obliquity was about 23°. 8, things were 
not very much different, and the same phenomenon took 
place about 22 days before and after winter solstice. Hav- 
ing the synthetic horizon, one can now ask what would be 
the azimuth of the setting sun on the solstice. This turns out 
to be 216°. 8 (disk center), and the horizon altitude is 1°.3, 
but in this case no reappearance is possible. As minimum 
and maximum azimuths for which the sun can shine on 
the back of the chamber are approximately 217° and 223° 



(Reijs 1998 1, both the solstice setting and the sun reappea- 
rance should be observable from within the chamber itself. 

The procedure outlined in this paper could be used, for 
instance, to test whether the reappearing sun is visible from 
other sites of archaeological relevance close to Maes Howe. 



* This correction is valid only on average, as the exact refraction de- 
pends on the physical conditions along the line of sight, and it is highly 
variable. See for instance Sampson et al. 120031. 
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The method can be applied during the exploratory phase 
of an archaeoastronomical site's study, either to plan direct 
on-site surveys or to test working hypotheses. When accu- 
racies of a few tenths of a degree are sufficient (as is cer- 
tainly the case for solar alignments), the results provided by 
the method can be used directly for the orientation analysis, 
saving a significant amount of time during the field work. 
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Fig. 5 Synthetic horizon computed for the Ring of 
Brodgar and apparent sun path for 6q=— 21°. 7 (top), 
-23°. 5 (middle), and -20°. 1 (bottom). 

As an example I analyze the case of the Ring of Brodgar] 
(3°. 2280 W, 59°. 0014 N). This stone circle, located between 
the Lochs of Stenness and Harray, has a diameter of about 
104 m, and is the third largest in the British Isles (Rug- 
gles |1999] l. As the calculations show, around ±20 days from 
winter solstice the sun sets over the Kame of Hoy, about 10° 
S of Ward Hill, without re-emerging from behind it (Fig. l5] 
top panel). At present day's winter solstice, the sun still sets 
behind the Kame of Hoy, without re-appearing (Fig.lSl mid- 
dle panel). Even for the increased obliquity expected for 
2700 BCE the apparent sun path does intersect the pro- 
file of Ward Hill. On the contrary, the reappearance of the 
sun is in principle observable from the Ring of Brodgar 
when i50=— 20°.l (November 22/January 21), that is about 
one month before and after winter solstice (Fig. |5] bottom 
panel). 

8 Conclusions 



The procedure I have outlined in this paper allows one to 
easily compute a synthetic horizon profile for any site. The 
accuracy in the final product depends on the resolution of 
the DEM and its vertical accuracy. Rms deviations of ^0°.l 
degrees are achievable with the SRTM90 data when the ho- 
rizon is at distances of the order of 10 km. For shorter dis- 
tances, the deviations are expected to be larger, especially if 
the horizon is defined by steep mountain profiles. In those 
cases the morphology of the area surrounding the site of in- 
terest needs to be examined more closely, and the synthetic 
horizon validated with direct theodolite measurements ob- 
tained along some critical directions. 



See http: //www. orkneyjar. com/hi story /brodgar/ 
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A Derivation of Earth's curvature 
correction 




In the following I will assume Earth can be described by an ellip- 
soid with semi-axis a and b. However, for the sake of simplicity, 
I will also assume that locally it can be approximated by a sphere 
having the curvature radius of the ellipsoid at the site under exam- 
ination. If (j) is the site geodetic latitude, this is given by 
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Throughout this paper I adopt the WGS84 ellipsoid, which has 
a=6378.1 km and flattening /= 1.0/298.25722356 (b=6356.8 km). 
If we now imagine to have an observer placed in A, at elevation 
/lA above the ellipsoid, looking at a site B placed at elevation Iib 
and with an angular separation on the ellipsoid defined as 9 (see 
Figure [ATJ, the altitude q of B as seen from A above the local 
horizontal plane is given by 



a = arctan — 
d" 
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While d" is the 



where d" and h" are defined as in Figure 
distance between A and B projected on the local horizontal plane, 
h" is the elevation of B above the same plane. The whole problem 
reduces to determining these two lengths. 

Looking at Figure [AT] we immediately note that: 

d" = {R + hs) sine . 



If we define A/i = HH' (see Figure Al i, then we have that: 

h — hs cos 8 — hA — Ah, . 

Since Ah — R{1 — cos 6), we finally obtain: 

h — hB cos9 — hA ~ R{i ~ cos9). 

Therefore, the elevation of B as seen from A is given by the 
following expression: 



a — arctan 
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(A3) 



This formula can be approximated considering that in the prac- 
tical applications 9 is going to be at most a few degrees. In fact, if 
we consider two points separated by an ellipsoidal distance ii=300 



Fig. Al Geometry of the problem. 



km, the angular separation is ^ ~ d/R=2° .1. In these circum- 
stances, one can use the following approximate expressions for 



the two trigonometric functions: cos9 ~ 1 — ^9 and sm9 ~ 9 
(where 9 is expressed in radians). 

After substituting them in Equation |A3| and considering that 
in all practical cases Hb/R ^1, we finally arrive at the following 
expression: 
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which can be readily used to derive the curvature-corrected alti- 
tude. 
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